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ABSTRACT 


^This  technical  report  describes  an  application  of  the  maximum  likelihood 
approach  for  the  identification  of  aircraft  parameters.  The  theory  and  the 
computational  aspects  of  the  likelihood  approach  are  discussed  in  detail. 


A computer  program  is  written  to  implement  this  techique.  The  details 
of  the  con^uter  program  are  given  in  a separate  document. 
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I.  INTRODUCTION 

To  meet  current  aircraft  and  rotorcraft  flight  testing  requirements, 
there  is  a genuine  need  for  an  accurate  procedure  to  evaluate  stability  and 
control  coefficients  from  recorded  response  data.  Most  techniques  developed 
in  the  past  have  either  not  had  the  capability  to  treat  highly  noisy  data  or 
could  not  be  practically  implemented  within  a con^juter  program  which  is  able 
to  identify  several  aircraft  parameters  at  a time.  The  work  reported  here 
and  the  resulting  conputer  program,  SCIDNT,  is  a development  to  fulfill  this 
requiranent. 

Past  methods  have  been  mostly  based  on  least  square  techniques  (e.g., 
curve  fitting  or  the  equation  error)  or  by  approximating  aircraft  motions 
by  low  order  systems  to  determine  the  frequency  and  damping  of  specific 
aircraft  modes.  The  resultant  algorithms  required  flight  test  time  that  may 
be  excessive  and  only  a few  parameters  can  be  identified.  The  basis  of  the 
present  approach,  on  the  other  hand,  is  the  likelihood  principle  of  statistical 
inference  theory.  The  likelihood  function  is  defined  and  the  parameter 
values  are  chosen  to  maximize  this  function.  This  approach  leads  to  an 
algorithm  which  resolves  many  of  the  problems  encountered  by  less  conpre- 
hensive  methods. 

This  report  is  a companion  document  to  the  User's  Guide  for  SCIDNT  [4J  , 
and  provides  the  SCIDNT  user  with  the  theoretical  foundations  of  the  program. 

In  particular,  the  optimization  methods  used,  the  Kalman  filter  used  for 
identification  in  the  presence  of  process  noise,  and  the  complete  sensitivity 
equations  are  discussed.  Examples  of  SCIDNT  performance  are  also  given. 

The  report  is  organized  as  follows.  Chapter  II  discusses  the  maximum 
likelihood  method  and  the  numerical  procedure  which  could  be  used  to  identify 
the  aircraft  stability  and  control  coefficients  from  flight  data.  Chapter  III 
presents  the  identification  results  from  simulation  data  as  well  as  flight 
data.  Here,  some  of  the  data  preprocessing  algorithms  which  were  used  to 
reduce  the  flight  data  are  also  given.  Chapter  IV  summarizes  the  report.  An 
expression  for  the  likelihood  function  is  derived  in  Appendix  A and  Appendix  B 
gives  the  aircraft  equations  of  motion. 
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II.  MAXIMUM  LIKELIHOOD  ESTIMATION  OF  AIRCRAFT  RARAMETERS 


2.1  imODUCTION 

The  maximum  likelihood  method  is  based  on  a fundamentally  different 
statistical  concept  from  least  square  methods  for  the  identification  of 
parameters  from  input- output  data.  Suppose  it  is  possible  to  make  a set 
of  observations  on  a systan,  vdiose  model  has  p unknovm  parameters  9.  For 
any  given  set  of  values  of  the  parameters  6 fron  the  feasible  set  0,  we 
can  assign  a probability  p(zl9)  to  each  outcome  z.  If  the  outcome  of  an 
actiial  experiment  is  z,  it  is  of  importance  to  know  which  sets  of  values  of 
9 might  have  led  to  these  observations.  This  concept  is  embedded  in  the 
likelihood  function  This  function  is  of  fundamental  importance 

in  estimation  theory  because  of  the  likelihood  principle  of  Fisher  and  others 
which  states  that  if  the  system  model  is  correct,  all  information  about  unknown 
parameters  is  contained  in  the  likelihood  function.  The  maximum  likelihood 
method  finds  a set  of  parameters  9 to  maximize  this  likelihood  function 

9 = max  ^(9/z)  C2.1) 

9e0 

In  other  words,  the  probability  of  the  outcome  z is  higher  with  parameters 
in  the  model  than  with  any  other  values  of  parameters  from  the  feasible  set. 
Conceptiially,  this  technique  can  be  summarized  as  follows: 

"Find  the  probability  density  functions  of  the  observations  for  all 
possible  combinations  of  unknown  parameter  values.  Select  the  den- 
sity function  whose  value  is  highest  among  all  density  functions  at 
the  measured  values  of  the  observations.  The  corresponding 
parameter  values  are  the  maximum  likelihood  estimates." 

Suppose  9 can  take  three  possible  values:  9^,  92 , and  9^.  Let  the 
probability  density  functions  of  observations  z for  these  three  values  of 
9 be  shown  in  Figure  2.1.  Then,  if  the  actual  observation  is  z,  92  is  the 
maximum  likelihood  estimate  of  9. 
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Figure  2.1  Maximum  Likelihood  Estimates 
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Usually,  it  is  more  convenient  to  work  with  the  negative  of  the  logarithm 
of  the  likelihood  function  (it  is  possible  to  do  so  because  the  logarithm  is  a 
monotonic  function).  Starting  from  a priori  values,  the  parameters  are 
updated  so  that  the  value  of  negative  log -likelihood  function  J at  the 
observed  value  of  the  outputs  decreases  continuously. 

The  great  asset  of  the  maximum  likelihood  method  is  that  it  can  be  used 
with  linear  or  nonlinear  models  in  the  presence  of  process  and  measurement 
noise.  The  maximum  likelihood  estimates  are  asymptotically  unbiased, 
consistent  and  efficient. 

The  details  of  the  maximum  likelihood  identification  procedure  are  given 
in  Appendix  A.  A simplified  flow  chart  is  shown  in  Figure  2.2.  Since  the 
Kalman  filter  is  a maximum  likelihood  approach  to  system  state  estimation, 
it  is  used  to  determine  the  state  time  history  with  fixed  parameters.  The 
negative  log -likelihood  function  which  depends  on  the  innovations  (difference 
between  the  observations  and  estimated  value  of  the  observations)  and  the 
covariance  of  the  innovations  can  then  be  determined.  The  first  and  second 
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^adients  of  J require  the  confutation  of  gradients  of  innovations  and  its 
covariance  with  respect  to  various  unknown  parameters.  These,  in  turn,  are 
determined  by  propagating  state  sensitivity  equations.  The  parameter  step 
size  which  reduces  J is  computed  using  modified  Newton-Raphson  or  a Gauss 
Newton  procedure.  Notice  that  it  is  possible  to  consider  noise  statistics 
as  parameters  and  identify  them.  A byproduct  of  this  numerical  procedure  is 
the  information  matrix.  The  information  matrix  provides  measures  of  the  identi- 
fiability  of  the  parameters.  If  this  information  matrix  is  singular  for  a 
certain  excitation  of  the  aircraft,  then  not  all  parameters  in  the  model  are 
identifiable  from  this  maneuver.  The  eigenvalues  and  eigenvectors  of  this 
symmetric  matrix  give  useful  information  about  how  well  the  parameters  can  be 
identified  from  given  data. 


2.2  MAXIMUM  LIKELIHOOD  IDENTIFICA. xON 

There  are  two  inf  lamentations  of  the  likelihood  method  which  have  been 
used  here. 

(a)  Output  error,  in  vdiich  the  noise  in  inputs  and  gust  effects  are 
small  and  can  be  neglected. 

(b)  Process  noise,  in  which  the  input  noise  and/or  gust  effects  are 
significant  and  must  be  included.  The  gust  and  input  noise  R'IS 
values  can  be  identified,  if  required. 

In  the  output  error  implementation,  the  state  equations  are  propagated 
directly  to  determine  the  estimated  state  vector  at  any  time.  The  Kalman 
filter  is  not  required  because  there  are  no  unknown  disturbances  in  the  state 
equations  which  would  cause  error  in  the  estimated  value  of  the  state. 
Therefore,  the  innovations  and  gradients  of  the  innovations  are  determined 
directly  by  integrating  the  state  equations  and  sensitivity  equations.  In 
the  process  noise  implementation,  the  state  equations  are  forced  by  unknown 
random  inputs.  A Kalman  filter  must,  then,  be  used  to  determine  the  estimated 
state  variable  vector.  The  gradients  of  the  innovations  for  system  parameters 
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are  computed  by  integrating  the  sensitivities  of  the  estimated  state  in  the 
Kalman  filter  equations.  The  process  noise  covariance  can  also  be  identified 
because  the  Kalman  filter  gains  and  innovations  covariance  depend  on  this 
parameter. 


The  maximum  likelihood  method  is  very  general  and  many  features  can  be 
added  to  any  implementation.  The  following  features  are  considered 
important  for  applications. 


(i)  Detailed  instrument  models  are  required.  The  measurements  are 
assumed  to  have  bias  and  random  noise  and  correction  terms  for 
coupling  with  other  states.  In  addition,  the  angle-of-attack 
and  sideslip  measurements  usually  have  scale  factor  errors. 

For  example,  the  angle-of-attack  measurement  is 

^m  ” q + ba  + no  (2.i;) 

where  ka  is  the  scale  factor,  ba  is  the  bias  and  no  is  random  noise. 


(ii)  In  most  implementations,  the  measurement  noise  covariance  matrix 
is  estimated  by 


R = 


N 


N 

.1 

i=l 


vCi)  v‘Ci) 


(2.2) 


where  v(i)  are  innovations.  In  general,  for  finite  data  length, 
this  would  give  an  estimate  of  measurement  noise  covariance  matrix 
with  all  nonzero  elements  and  no  special  structure.  In  other 
words,  Eq.  (2.2)  assumes  that  there  are  m(m+l)/2  unknown  elements 
in  measurement  noise  statistics.  If  the  structure  of  R is  known, 
as  is  usually  the  case  in  aircraft  applications,  the  number  of 
unknown  eloments  in  R is  fewer  than  m(m+l)/2.  For  example,  if 
the  measuranent  noise  in  different  channels  is  independent,  only 
diagonal  terms  of  R need  to  be  estimated.  Using  Eq.  (2.2)  will 
lead  to  an  overparameterized  model  and  incorrect  parameter  estimates. 
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The  statistically  correct  estimate  of  measurement  noise  covariance 
matrix  for  independent  measurement  noise  sources  is, 

1 ( N ™ 

R = diag  ^ v(i)  v (i) 


(2.3) 
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Information  Ag.gregation 

.Any  available  a priori  information  about  parameters  can  be  incorporated 
into  the  identification  program.  This  is  done  by  finding  the  weighted  mean 
of  the  a priori  parameter  values  and  the  output  of  the  identification  program, 
the  weights  varying  inversely  as  the  corrresponding  covariances.  This  has 
been  done  before  [1] . A similar  approach  can  be  used  to  combine  estimates 
from  different  identification  runs  and  thus  make  full  use  of  available  data. 

In  the  vector  case,  let  9^  and  02  be  two  different  estimates  of  the 
parameter  vector  with  information  matrices  and  M2  respectively.  Then  we 
can  aggregate  these  two  estimates  to  get  a new  information  matrix  and 
parameter  vector. 


M = M^  + M2 

(2.5) 

§ = M ^ ^2^2^ 

(2.6) 

2.3  ODMPUTATIONAL  ASPECTS  OF  MAXIMUM  LIKELIHOOD  ESTIMATES  [2] 

The  application  of  maximum  likelihood  identification  in  practice  requires 
efficient  computational  algorithms  to  maximize  the  likelihood  function.  It 
is  not  uncommon  to  find  situations  where  the  likelihood  surface  has  multiple 
maxima,  saddle-points,  discontinuities  and  singular  Hessian  in  the  parameter 
space.  The  application  of  the  steepest  descent  method  leads  to  extremely 
slow  convergence  rate  and  the  straightforward  application  of  the  Newton- 
Raphson  and  the  Gauss -Newton  methods  may  lead  either  to  no  convergence  or 
convergence  to  wrong  stationary  points.  From  a statistical  veiwpoint,  only 
the  maximum  of  the  likelihood  function  provides  an  unbiased,  consistent  and 
efficient  estimate.  Thus,  it  is  important  to  locate  the  true  maximum  of  the 
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likelihood  function.  With  the  present  computation  methods  of  nonlinear 
programming,  in  general,  this  could  be  an  extremely  difficult  and  time-consuming 
task.  In  aircraft  applications,  a priori  information  about  the  parameters 
is  used  to  choose  a good  starting  value  and  impose  physically  meaningful 
constraints  on  the  parameters.  This  increases  the  chances  of  approaching 
the  absolute  maxima. 

2.3.1  Identifiability  Problems  and  Solutions 

Anomalies  also  occur  in  the  likelihood  function  due  to  inadequate  model 
specification  and  parameterization.  Sane  of  these  anomalies  that  lead  to 
numerical  difficulties  are:’ 

1.  Discontinuities  and  Singularities:  There  are  examples  where  the 
likelihood  function  increases  up  to  a maximum  value  and  then  falls 
inmediately  to  minus  infinity.  In  other  cases,  the  likelihood 
function  may  rise  suddenly  to  infinity  for  some  particular  parameter 
values.  These  problems  can  often  be  avoided  either  by  changing  the 
model  or  by  changing  the  parameterization. 

2.  Singular  or  Nearly  Singular  Information  Matrix:  This  is  one  of  the 

most  common  problems  in  parameter  identification  since  it  is  difficult 
to  determine  a priori  whether  a given  sample  has  adequate  information 
for  estimating  all  the  parameters  in  the  model.  Generally,  the 
tendency  in  practical  applications  is  to  include  all  the  parameters 
about  which  there  is  some  uncertainty.  This  leads  to  overpara- 
meterization which  in  turn  produces  a singular  or  nearly  singular 
information  matrix.  The  likelihood  surface  in  turn  has  long, 

curved  and  narrow  ridges  along  which  the  convergence  rate  is 
extremely  slow.  Under-parameterization  does  not  really  provide 
the  solution  since  it  may  lead  to  spurious  local  maxima  and 
saddle-points. 
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The  basic  iteration  in  gradient -type  nonlinear  programming  methods  is 


9.  - = 0.  - p.  RT-^g. 

X+1  1 1 1 


(2.7) 


.th 


where  9.  is  the  parameter  vector  at  the  i iteration,  g.  is  a vector  of 
gradients  of  the  negative  log- likelihood  function  J(9),  i.e.,  the  j 
component  of  g^  is  given  as 

ssfor)!  c.i 


(2.8) 


is  an  approximation  to  the  second  partial  matrix 


where  is  a scalar  step  size  parameter  chosen  to  ensure  that  J(9^^j^)  < J(0^)  - e 
where  e is  a positive  number  that  can  be  chosen  in  a variety  of  ways.  In 
the  Gauss-Newton  method  used  here,  the  matrix  of  second  partials  is  approximated 
by 


RlO.k)  = M.Cj.k)  ^ 1(5*^  8-1 


. 1.  n._  0-1  36  0-1  36  I - 

^ Tr  B B \ 


(2.9) 


where  M is  the  approximation  of  R.  The  first  partials  of  innovations  v and 
their  covariance  can  be  computed  by  solving  linear  difference  equations,  i.e., 
the  sensitivity  equations.  Notice  that  in  this  approximation  computation  of 
second  partials  of  state  and  covariance  is  not  required. 
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There  are  several  problems  with  this  method.  As  long  as  is  positive 
semi-definite,  it  is  possible  to  find  a step  vdiich  will  reduce  the  cost. 

As  mentioned  before,  most  of  the  problems  arise  vdien  the  information  matrix 
is  almost  singular,  i.e.,  the  eigenvalues  are  spread  far  apart.  This  near- 
singularity  of  the  information  matrix  manifests  itself  in  several  ways. 

(a)  rXiring  the  computation  of  the  inverse  of  the  information  matrix, 
the  positive-definiteness  may  be  lost  because  of  round-off  errors. 


(b)  The  step-size  may  be  very  large  in  those  directions  which  correspond 
to  small  eigenvalues  of  the  information  matrix.  Let 


in  m 

M = y X.V.V. 
i=l  ^ ^ " 


(2.10) 


m 


(2.11) 


^ 1 T 

IT  VlS 


■ - v-1 


(2.12) 


where  are  the  eigenvalues  and  v.  the  corresponding  eigenvectors  of  M.  Thus, 

X X 

if  v|g  is  not  small  for  small  eigenvalues  X.,  the  step  in  the  direction  of  v. 
is  large.  Usually,  for  small  eigenvalues,  vig  is  also  small,  but  it  is  the 
difference  of  large  and  almost  equal  numbers.  Therefore,  in  most  cases,  it  is 
only  a numerical  problem.  Coiqjlications  arise  from  local  irregularities  in 
the  likelihood  function.  Several  remedies  can  be  used,  depending  on  the 
problem  complexity. 

(a)  A9  * -oM 


Let  A6'  = M'^g 


(2.13) 


then  A0 


-aA0' 


(2.14) 
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Positive-definiteness  of  M could  be  lost  during  the  inversion. 
This  can  be  avoided  by  finding  A0'  by  solving  the  following 
linear  equations 


MA9"  = g 


(2.15) 


The  accuracy  of  Eq.  (2. IS)  also  depends  on  conditioning  of  M 
(ratio  of  maximum  and  minimum  eigenvalues)  but  will  be  much 
superior  as  compared  to  Eq.  (2.13). 


(b)  Conditioning  of  M can  be  improved.  Premultiply  both  sides  of 
Eq.  (2.15)  by  A where 


A - diag 
to  get 

(A'Hia'^)  (AA0')  = A*^g 


(2.16) 


or 


M*A0*  = g* 


(2.17) 


All  elements  of  M*  are  between  -1  and  +1.  In  general,  this  matrix 
is  much  better  conditioned  than  M.  Equation  (2.17)  can  be  used 
to  solve  for  A0*  and  then  A9'  is  obtained  by  rescaling. 


(c)  The  rank  deficient  inverse  of  M can  also  be  used.  In  this 

procedure,  the  eigenvalues  of  matrix  M are  arranged  in  decreasing 
order  of  magnitude  and  all  eigenvalues  below  a threshold  b are 
neglected  in  using  Eq.  (2.11)  for  the  inverse,  i.e. , 


(2.18) 
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where 


1 

I 


i 


r 


X.  > b i < k 

1 — — 


X^  < b i > k (2.19) 

This  is  equivalent  to  searching  for  the  maximum  of  the  likelihood 
function  in  the  subspace  spanned  by  v^(i  £ k) . If  this  process 
is  continued  until  the  end,  the  f'Stimates  may  never  converge  to 
their  maximum  likelihood  values.  It  should  also  be  noted  that 
since  only  large  eigenvalues  are  included,  step  A6  may  be  "small". 

(d)  Instead  of  neglecting  eigenvalues  of  M,  vdiich  are  smaller  than  a 
certain  threshold  b,  these  smaller  eigenvalues  are  increased  to  b. 
Thus,  the  inverse  of  M is  approximated  as 


k ^ rp  T ni 
y 1 v.vT  . 1 S 

1-1  H ^ i4i 


T 

v.v. 
1 1 


(2.20) 


In  this  case,  when  the  program  converges,  the  likelihood  function 
will  reach  a maxima.  Convergence  may  be  slow. 


(e)  This  fix  of  the  Gauss-Newton  method  is  a two  step  procedure.  The 
step  is  taken  by  finding  M ^ using  Eq.  (2.18).  Then,  another  step 
direction  is  computed  using  only  the  small  eigenvalues,  i.e., 

T 

m v.v.g 

A9,»-a  1 (2.21) 

^ i=k+l  ^i 


Then  step  size  a is  chosen  by  a one  dimensional  search.  As  the 
algorithm  approaches  the  maximum,  it  may  be  worthwhile  to  carry 
out  one  dimensional  search  in  each  v^(i  > k)  direction. 

The  last  approach  when  used  with  normalized  information  matrix  NI* 
of  Eq.  (2.17)  is  one  of  the  best.  It  is  recommended  if  the  number 
of  unknown  parameters  is  large.  In  most  cases,  fix  (d)  is  adequate. 
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III.  EXAMPLES 
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3.1  INTRODUCTION 


The  computer  program  (SCIDNT)  is  based  on  the  maximum  likelihood 
method  outlined  in  Chapter  II  and  Appendix  A.  Many  examples  have  been 
performed  to  test  the  usefulness  and  accuracy  of  the  technique.  These 
exan^iles  include  both  longitudinal  and  lateral  motions  and  simulation  and 
flight  test  data  from  a high  performance  swept-wing  Navy  fighter. 


The  equations  of  motion  used  in  the  identification  are  given  in  Appendix 
B.  In  addition  to  a,  u,  q and  9 there  is  an  additional  state  in  the  longi- 
tudinal case  corresponding  to  the  gust  angle-of-attack,  Similarly, 
there  is  a fifth  state  in  the  lateral  equations  of  motion  which  corresponds 
to  random  fluctuations  in  sideslip  angle  because  of  lateral  gusts. 

3.2  SIMULATION  DATA 

Simulated  data  provides  a tool  for  program  checkout  and  experimentation 
for  special  problems  (e.g. , input  design).  In  this  section,  the  important 
simulated  data  runs  by  SCIDNT  are  summarized.  A DC-8  simulation  at  600  fps 
forward  speed  and  sea  level  altitude  was  used  both  for  longitudinal  and 
lateral  motions.  The  computer  program  PLANE  generates  the  simulation  data. 

3.2.1  Longitudinal  Motion 

The  values  of  parameters,  used  in  the  simulation  to  generate  the  data,  are 
shown  in  Table  3.1  An  elevator  input  is  used  to  excite  the  aircraft  [see 
Figure  3.1].  Five  noisy  instnjments  measure  the  resulting  aircraft  response. 

It  is  clear  that  only  the  short  period  mode  is  excited. 

In  the  first  case,  it  is  assumed  that  there  are  no  gust  effects.  SCIDNT 
is  started  using  least  square  parameter  values.  Five  parameters,  Za,  Nb,  Mq, 

and  Mjg,  vrfiich  determine  the  short  period  motions  are  identified  in 
addition  to  bias  and  random  noise  in  instruments.  The  identified  parameters 
are  shown  in  Table  3.1.  Also  shown  are  the  root  mean  square  values  of  the 
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Table  3.1 


X-8  SIMULATED  DATA 

Velocity  600  fps 
Altitude  Sea  Level 


Parameter 

True  Value 

Estimated  P 
fvdth  1 

arameter  Value 
a bound] 

Output  Error 

Process  Noise 

Z 

a 

-0.3237 

-0.3168 

-0.3205  (.0025) 

Z 

-1.074 

* 

* 

u 

Z 

q 

1.0 

0.9212 

1.0313  (.057) 

^0 

0.0 

* 

* 

X 

0.02845 

* 

* 

a 

-0.006386 

* 

* 

X 

0.0 

* 

* 

q 

M 

-1.204 

A 

* 

a 

M 

0.0 

A 

* 

u 

M 

q 

-0.3735 

-0.5006 

-0.402  (.012) 

Z6e 

-0.0178 

-0.131 

0.0  (.056) 

^6e 

0.0 

* 

* 

^6e 

-1.381 

-1.411 

-1.396  (.014) 

0) 

0.5 

* 

* 

c 

Q 

0.002 

* 

0.0024  (.00018) 

k 

1.0 

* 

* 

a 

ka£a 

0.0 

* 

fc 

V 

0.005 

* 

0.011 

a 

0.5 

* 

0.93 

u 

0.005 

A 

0.0105 

n 

0.005 

* 

0.0102 

0.5 

* 

0.5 

* Not  identified 

[All  in  units  of  ft.,  rad.,  sec.] 
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Figure  3.1  Elevator  Input  for  DC-8  Simulation 


random  noise  in  instruments.  The  time  history  of  the  simulated  and  predicted 
aircraft  response  are  shown  in  Figure  3.2.  The  fit  is  good.  Next,  a simulation 
data  is  generated  with  gust  effects  included.  The  pwoer  spectral  density  of 
the  white  noise  driving  the  gust  state  is  .002.  Five  short  period  parameters 
and  the  power  spectral  density  of  the  white  noise  are  estimated  from  the 
simulation  data.  The  time  history  fits  are  given  in  Figure  3.3  and  the 
estimated  values  in  Table  3.1.  Also  shown  are  one  standard  deviation  values 
on  parameter  estimates. 


t 

* 

? * 
► ^ 


It  is  in^Jortant  to  understand  the  significance  of  the  standard  deviation 
on  parameter  estimation  errors.  Because  of  instrument  errors  and  unknown 
gust  effects,  the  estimated  parameter  values  are  not  the  same  as  true 
parameter  values.  In  other  words,  there  is  an  estimation  error.  Consider 
the  parameter  estimates  to  be  random  samples  from  a normally  distributed 
population  whose  mean  (unknown)  is  the  true  parameter  value.  Then  the  parameter 
estimates  are  within  two  standard  deviations  of  the  true  parameter  value 
about  95%  times.  Usually,  the  parameter  estimates  do  not  have  a normal 
distribution  and  a relation  such  as  the  Chebyshev  inequality  must  be  used  to 
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give  more  useful  meaning  to  the  standard  deviation*.  Notice  that  this  is  a 
probabilistic  concept  and  does  not  ensure  anything  definite  in  a specific  case. 
The  maximum  likelihood  method,  when  used  with  the  numerical  procedure  of 
Section  2,  gives  the  standard  deviations  of  parameter  estimation  errors  as  a 
byproduct  of  the  parameter  identification. 

To  explain  the  meaning  of  the  estimation  error  standard  deviation  in  this 
example,  consider  the  parameter  2a.  Assume  that  we  do  not  know  the  true  value 
of  this  parameter  used  in  the  simulation.  In  the  process  noise  case  the 
estimated  value  is  -.3205.  We  know  it  is  not  the  true  Za,  but  is  some  value 
close  to  the  true  value.  What  additional  information  does  a standard  deviation 
of  .0025  give  us?  If  the  Chebyshev  inequality  is  used,  we  can  be  75%  certain 
that  the  true  parameter  lies  within  +.3205  + .0050,  i.e.,  -.3155  and  -.3255. 

3.2.2  Lateral  Motions 


A rudder  input  followed  by  an  aileron  input  of  Figure  3.4  are  used  to 

excite  the  lateral  motions  of  DC-8  aircraft.  The  parameter  values  used  in 

the  simulation  are  shown  in  Table  3.2.  Only  five  parameters,  N , N^,  Yq, 

r i3  p 

and  N5  are  identified  while  the  others  are  fixed  at  the  least  square 
values.  The  identified  value  and  the  least  square  starting  value  for  the 
maximum  likelihood  program  are  also  shown  in  Table  3.2. 

The  estimated  value  of  Ng,  Yg,  and  is  quite  close  to  the  true 

value.  However,  there  is  an  error  of  18.7%  in  the  identified  value  of  N . 

v 

This  does  not  necessarily  mean  that  is  poorly  identifiable.  The  error 
could  occur  because  other  parameters  have  been  fixed  at  incorrect  values 
in  the  maximum  likelihood  identification. 


* The  Chebyshev  inequality  states  that  a random  variable  will  lie  within 
+XO  from  the  mean  more  than  (1  - 1 ) fraction  of  times.  For  e.xample,  a 

? 

random  variable  will  lie  within  ^2a  more  than  75%  times.  Note  that  this 
inequality  is  independent  of  the  distribution  function. 


Table  3.2 


RARAMETER  IDENTIFICATION  RESULTS  FOR  LATERAL  CASE  EXA^IPLE 

^.(SIMULATED  DATA) 


RARANETER 

TRUE  VALUE 

STARTING 

VALUE 

rVA  RESULT 

ML  RESULT 

% ERROR 

L 

n 

-2.029 

-12.840 

-2.105 

* 

L 

r 

0.4697 

- 2.700 

.4407 

* 

L 

-2.092 

1.628 

-1.883 

* 

N 

n 

-0.05828 

- 7.770 

-0.05089 

* 

N 

r 

-0.2724 

-17.540 

-0.3107 

-0.3233 

13. ' 

^3 

1.655 

1.360 

1.716 

1.693 

2.3 

^3 

-0.1274 

-20.1896 

-0.1289 

-0.1348 

5.7 

a 

0.0 

1.186 

0.715 

0.005 

^6a 

22.02 

- 1.83 

22.39 

21.98 

0.2 

0.4807 

- 0.260 

0.2576 

* 

-0.1709 

1.0534 

-0.2800 

* 

N^r 

-1.3680 

- 0.0071 

-1.3734 

-1.3116 

4.1 

^5a 

0.0 

0.2314 

-1.274 

* 

Y6r 

0.03630 

0.0114 

0.23232 

it 

* Not  identified 

[All  in  units  of  ft.,  rad.,  sec.] 


The  predicted  measurement  time  histories  based  on  identified  parameters 
are  a good  fit  to  the  "true"  measurements,  as  shown  in  Figure  3.5. 

3.3  FLIGHT  TEST  DATA 


The  flight  test  data  presented  in  this  section  is  from  a Navy  swept-wing 
fighter  at  30,000  feet  altitude,  0.6  mach  airspeed  and  20°  sweep  angle.  The 
wind  tunnel  estimates  of  aircraft  stability  and  control  coefficients  are 
inaccurate. 
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ExnnPlE  LStERBL  CASE  MITM  nERSUREnENT  NOISE  tXBnPLE  lHIEKAL  CRSE  MUH  hEBSUREBENI  NOISE 


3.3.1  Longitudinal  Motions  Flight  Data 


The  response  data  shown  in  Figure  3.7  is  for  a doublet  input  of  Figure  3.6. 
The  data  is  fairly  noise  free.  "Spikes",  arising  from  instrument  recording 
dynamic  response  characteristics  are  evident. 

Since  no  a priori  parameter  estimate  of  sufficient  accuracy  to  start  the 
maximum  likelihood  identification  program  were  available,  it  was  necessary  to 
use  a least  squares  approach  to  obtain  starting  values  of  the  parameters. 

These  values  are  shown  in  Table  3.3.  The  flight  test  data  is  processed  as 
such  to  obtain  the  parameter  estimates  of  Table  3.3  and  predicted  measurement 
time  history  of  Figure  3.7.  The  estimates  look  reasonable  and  the  fit  to  the 
actual  measurement  is  quite  good.  The  outliers  which  are  caused  by  recording 
or  instrumentation  error  can  cause  significant  errors  in  parameter  estimates. 
Therefore,  the  data  is  preprocessed  before  identification  is  performed.  In 
this  preprocessing  the  outliers  are  ranoved.  In  addition,  one  more  parameter 
considered  important  is  identified.  The  identified  parameter  values  with  data 
preprocessing  is  given  in  Table  3.3  and  time  history  fits  in  Figure  3.8. 

There  is  a small  change  in  most  estimates.  The  one  standard  deviation  esti- 
mation errors  are  enclosed  in  the  brackets. 


FM  OflTfl. OUTPUT  ERROR 


Figure  3.6  Elevator  Input  for  Flight  Test  Data 
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Figure  3.7  (cont'd)  Flight  Test  Data  With  Output  Error  (Outliers  Not  Re?noved) 


Table  3.3 


F-14  RUN2 

Altitude  30,000  ft. 
Speed  362  knots 
Sweepback  20“ 


DERIVATIVE 

STARTING  VALUE 
(LE4ST  SQUARES 
.AND  WIND  TUNNEL) 

ESTIMATED 
RA!V  DATA 

VALUE  OVim  la  BOUND) 
OUTLIERS  REMOVED 

Z 

-0.38 

-0.347  (.033) 

-0.339  (.00271) 

Z 

-0.00066 

* 

A 

u 

Z 

1.04 

* 

A 

q 

Yo 

0.0 

* 

A 

X 

285.0 

A 

A 

a 

Xu 

-3.308 

* 

A 

42.34 

A 

A 

M 

a 

-0.755 

-1.276  (.0789) 

-1.285  (.00649) 

Mu 

0.00416 

A 

A 

-0.600 

-0.736  (.179) 

-0.7068  (.011) 

“c 

0.50 

A 

A 

2«e 

-0.20 

0.0378  (.218) 

0.119  (.020) 

X6e 

0.0 

A 

A 

-5.46 

-5.357  (.452) 

-5.338  (.033) 

Q 

0.0 

A 

A 

k 

1.0 

A 

A 

a 

koto 

V 

0.0 

A 

-0.0754  (.0062) 

ba 

0.0 

A 

0.000124 

bu 

0.0 

A 

-0.0934 

bq 

0.0 

A 

0.00056 

'’e 

0.0 

A 

0.000213 

baz 

0.0 

A 

-0.05 

0.0075 

0.0428 

0.00159 

100.0 

52.0 

4.72 

0.007 

0.0424 

0.0038 

0.015 

0.0877 

0.0021 

0, 

0.9 

8.1 

0.77 

at 

* Not  identified 


[All  in  units  of  ft.,  rad.,  sec.) 

[Data  Supplied  by  NATC  TO  SCI  - .^ril,  1974] 
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Figure  3.8  Flight  Test  Data  (Outliers  Raroved) 


r 


1 

I 

i 


I 

j 


Another  response  record  resulting  from  a doublet  input  at  approximately 
the  same  flight  condition  was  obtained  three  months  later.  The  identified 
value  from  previous  run  were  used  for  startup.  ^Vhen  the  data  was  processed 
as  such,  it  was  immediately  clear  that  the  pitch  rate  measurement  is  bad. 

The  raw  data  was  rerun  and  the  aircraft  stability  and  control  derivatives 
were  identified  without  using  the  q measurement,  but  using  the  measurements 
of  a,  9,  u and  nominal  acceleration.  The  results  are  given  in  Table  3.4. 
Subsequently,  the  measurement  of  q was  corrected  (it  had  several  sign  errors) 
and  also  included  in  the  identification  to  obtain  parameter  values  shown  in 
the  last  column  of  Table  3.4.  The  time  history  plots  of  the  true  and  predicted 
measurements  are  given  in  Figure  3.9.  Some  parameters,  e.g.,  Za  and  Ma,  are 
quite  different  from  the  last  case. 

3.3.2  Lateral  Motions  Flight  Data 

Lateral^response  data  was  obtained  for  rudder  and  aileron  inputs  of 
Figure  3.10.  “’^e  data  was  processed,  first,  using  the  least  squares  approach 
to  obtain  the  starting  values  for  maximum  likelihood  identification.  These 
values  are  used  in'  tl^  output  error  mode  of  the  maximum  likelihood  method. 

This  Results  in  parameter  estimates  of  Table  3.5  and  time  history  plots  of 
Figure  3.11.  Next,  the  same  data  is  processed,  assuming  lateral  gusts.  The 
reuslts  are  shown  in  Table  3.5  and  Figure  3.12.  The  power  spectral  density 
of  the  process  noise  is  .000347  sec  Since  the  gust  model  break  frequency 
is  .5  sec,  this  gives  a root  mean  square  sideslip  gust  of 


(8  JRMS  = 

o 


.000347 
2x  5 


.01863  rad 


(3.1) 


The  aircraft  speed  is  362  knots.  The  lateral  wind  has  a root  mean  square 
value  of  7.68  ft  sec 

3.3.3  Verification  of  Parameter  Estimates  Obtained  from  Flight  Data 

The  following  steps  are  used  to  verify  the  validity  of  parameter 
estimates  obtained  from  flight  data. 

J 
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Table  3.4 


RUNl  (DOUBLET  INPUT) 


Altitude  30,000  ft. 
Speed  360  knots 
Sweepback  20® 


t 


; 

V 


* 


PARAMETER  VALUE  (WITH  la  ERROR  BOUND  IN  BRACKETS) 

PARAMETER 

STARTING  VALUE 
(FROM  PREVIOUS  RUN) 

q FRai  DATA 
(INCORRECT) 

q NOT  USED 

q CORRECTED 
A.ND  USED 

-0.3394 

0.545 

-0.4268  (.0097) 

-0.4089  (.0059) 

Zu 

-0.00066 

* 

A 

1.04 

A 

A 

A 

^0 

0.0 

A 

A 

A 

285.0 

A 

A 

A 

X 

-3.3 

A 

A 

A 

U 

X 

42.3 

A 

A 

A 

q 

-1.28 

-4.59 

-2.002  (.018) 

-1.88  (.0078) 

”u 

0.0042 

A 

A 

A 

-0.74 

-2.73 

-0.4907  (.024) 

-0.6726  (.0089) 

“*c 

O.S 

A 

A 

A 

2^e 

x<s. 

0.119 

0.0 

0.639 

A 

0.0856  (.014) 

-0.0147  (.0136) 

M«e 

S.34 

0.416 

4.154  (.093) 

5.30  (.033) 

Q 

0.0 

A 

A 

A 

k 

1.0 

A 

A 

0.881  (.0089) 

a 

kata 

V 

0.0 

-0.08 

-0.0594  (.0066) 

A 

b 

0.0 

-0.0328 

0.00609 

-0.00071 

a 

b 

0.0 

-2.06 

0.344 

-0.147 

u 

bp 

0.0 

-0.0104 

0.00164 

-0.000877 

H 

‘’e 

0.0 

-0.0378 

0.0112 

-0.00505 

0.0 

-0.81 

0.922 

-0.528 

0.007S 

0.034 

0.0072 

0.00262 

a 

100.0 

2.7 

4.05 

1.35 

u 

Op 

0.007 

0.023 

0.054 

0.0034 

M 

0.015 

0.063 

0.0131 

0.0100 

0.9 

10.0 

2.6 

2.1 

• Not  identified 

(All  in  units  of  ft.,  rad.,  sec.) 

(Data  Supplied  by  NATC  TO  SCI  - July,  1974] 
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Table  3.5 


UTERAL  MOTIONS  (RUNS)* 


Altitude  30,000  ft. 
Speed  360  knots 
Sweepback  20“ 


t 


STARTING  VALUE 
(LEAST  SQUARES) 

MAXIMJM  LIKELIHOOD  ESTD-IATE  (IVITH  lo  BOUND) 

PARAMETER 

NO  PROCESS  NOISE 

?R(X:ESS  NOISE 

Lp 

-1.975 

-2.004  (.0235) 

-1.948  (.0205) 

Lr 

4.877 

2.686  (.0645) 

2.369  (.0541) 

S 

-8.29 

-8.500  (.0982) 

-6.952  (.0846) 

-0.004 

* 

A 

-0.62 

* 

A 

1.96 

A 

* 

-0.05 

-0.0649  (.00190) 

-0.0819  (.00121)  ! 

0.0 

H 

A 1 

0 

0.0 

* 

A 

0 

10.29 

n 

A 

0.24 

It 

A 

0.0 

it 

A 

-2.828 

-3.337  (.0865) 

-1.446  (.0661) 

^6r 

1.S5 

* 

A 

0.42 

-0.015  (.0018) 

-0.0150  (.00121) 

1.0 

it 

A 

0.03 

A 

A 

“c 

0.5 

A 

A 

Q 

0.0 

A 

0.000347  (.000026) 

bp 

- 

-0.000209 

-0.00264 

- 

-0.000126 

0.000318 

*’8 

- 

-0.0000590 

0.000636 

- 

-0.00230 

-0.00516 

bay 

- 

0.00568 

-0.0629 

0.005 

0.0405 

0.0246 

a 

0.005 

0.0196 

0.0069 

y 

0.005 

0.0190 

0.00487 

0.005 

0.0192 

0.0245 

0.5 

1.57 

0.85 

* Not  identified 

[All  in  units  of  ft.,  rad.,  sec.) 

(Data  Supplied  by  NATC  TO  SCI  - July,  1974] 
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Figure  3.12  (cont'd)  Measured  and  Predicted  Aircraft  Responses  in  the  Lateral  Direction  (With  Process  Noise) 


(a)  The  time  history  matches  between  the  actual  measurements  and 
measurements  predicted  from  estimated  parameter  values  are  good. 

(b)  The  parameters  are  used  to  predict  the  frequency  and  damping  of 
various  modes.  These  are  con^ared  with  values  read  from  the 
response  of  different  inputs.  See  Tables  3.6  and  3.7  for 
longitudinal  and  lateral  motions,  respectively. 

(c)  Response  data  from  different  inputs,  but  at  the  same  flight 
condition  is  processed.  The  identified  parameter  values  from 
different  runs  are  close  to  each  other;  see  Table  3.8  for 
the  lateral  case. 

(d)  Identified  parameter  value  from  response  to  one  input  are  used 
to  predict  response  to  another  input  at  the  same  flight 
condition.  Again,  excellent  results  were  obtained. 
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IV.  CONCLUSIONS 

-I 


This  report  outlined  the  theoretical  and  statistical  background  of  the 
maximum  likelihood  method  for  parameter  estimation  and  showed  its  usefulness 
in  identifying  aircraft  stability  and  control  coefficients  frcsn  flight 
response  data.  Several  important  conclusions  can  be  drawn  regarding  the 
application  of  maximum  likelihood  method  and,  in  particular,  the  computer 
program  SCIIWr  for  reducing  flight  test  data. 

It  is  possible  to  obtain  meaningful  parameter  estimates  from  well 
conducted  flight  tests  using  the  maximum  likelihood  approach.  These  estimates 
are  consistent  between  runs,  can  predict  frequency  and  damping  of  aircraft 
natural  modes  and  aircraft  response  to  similar,  but  not  same,  inputs  at  the 
same  flight  condition  (Section  3.3.3). 


Both  the  output  error  and  the  process  noise  options  of  the  maximum 
likelihood  are  important.  While  in  most  runs,  the  output  error  technique 
gives  accurate  parameter  estimates  and  time  history  matches,  in  some  runs 
significant  improvement  is  achieved  with  the  inclusion  of  process  noise  in 
the  model.  See  Sections  3.2.1  and  5.3.2. 


Because  spurious  signals  are  often  present  in  flight  data,  some  data 
preprocessing  will  usually  be  required.  It  is  best,  therefore,  to  make  a 
visual  examination  of  the  recorded  data  before  using  any  identification 
algorithm.  The  data  should  be  checked  for  sign  inconsistencies,  outliers 
and  other  obvious  errors.  The  dimensions  of  various  signals  could  be  checked 
for  compatibility.  The  data  preprocessor  could  vary  in  con^jlexity 
(Section  3.3). 


« 

Input  signals. 

♦ 

with  care.  The  set 

t 

• 

with  which  they  can 

>•  ^ . 
* ** 
f 

in  Table  3.8) . 

used  to  excite  the  aircraft  response,  should  be  chosen 
of  parameters  which  can  be  identified  and  the  accuracy 
be  identified  depends  on  the  input  (compare  runs  6 and  8 
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In  summary,  the  computer  implementation  SCIDNT  of  the  maximum  likelihood 
parameter  identification  approach  should  prove  useful  in  reducing  flight  test 
data  to  determine  aircraft  stability  and  control  derivatives  on  a routine 
basis.  It  is  important  to  understand,  however,  that  the  flight  test  must  be 
planned  with  the  objectives  in  mind  that  the  response  data  must  be  carefully 
collected  and  checked  out  for  successful  evaluation  of  meaningful  and  useful 
aircraft  stability  and  control  coefficients. 


I- 

t 


I 
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APPENDIX  A 


MAXIMUM  LIKELIHOOD  IDENTIFICATION  OF  PARAMETERS 
A.l  GENERAL  NONLINEAR  SYSTBIS 

Consider  the  general  nonlinear  aircraft  equations  of  motion 

X = f(x,u,0,t)  + r(0,t)w  0 < t < T 

E(x(0))  = x^(0) 

E{(x(0)  - x^(0))(x(O)  - X (0))’’^}  = P (0) 
o o o 

where 

x(t)  is  n X 1 state  vector 

u(t)  is  2 X 1 input  vector 

0 is  p X 1 vector  of  unknown  parameters 

r is  n X q process  noise  distribution  matrix 

w(t)  is  q X 1 random  process  noise  vector 


Sets  of  m measurements 


yCt^) 


are  taken  at  discrete  times  t. 


y(tj^)  = h(x(tj^),u(tj^),0,tj^)  + v(tj^) 


k = 1,2,3 N 


w(t)  and  vCt^)  are  Gaussian  random  noises  with  the  following  properties 


«• 

t 

« 

► ' . 


E(w(0)  = 0 


E(v(tj^))  = 0 E(w(t)  V(tj^))  = 0 


E(w(t)  w (r))  = Q(0,t)  6(t-t),  E(v(t.)  V^(t,  ))  = R(0,t.)6.. 

J K J 


CA.l) 


(A.  2) 


(A.  3) 


(A.  4) 
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The  unknown  parameters  are  supposed  to  occur  in  the  functions  £ and 
h and  in  matrices  P,  Q,  R,  and  x^.  In  the  analysis  to  follow,  the  model 
and  the  functional  form  of  f and  h is  assumed  known  correctly. 

The  set  of  observations  y(t^)  , y(t2)  , . . . ,y  (tj^)  constitutes  the  outcome 
Z in  this  case.  The  likelihood  function  for  0,  which  has  the  same  form  as  the 
probability  of  the  outcome  z for  a certain  value  of  parameters  0,  is  given  by 

i^(0|z)  = p(z|0) 


= p(y(tp.  y(t2)  — 

= P(Yj^|0)  (A. 5) 

where 


\ ~ ^ ^ 

p(Yj^|9)  = p(y(tN^lVi'®^*p(Vi!®^ 

= pCyC^lYj^.^.e)  pCy(tj^.i)|Yj^.2.e)  p(Yj^_2ie) 

N 

= _n  p(y(t.)|Y._j,0)  (A. 6) 

The  log -likelihood  function  is 

N 

log  (5?(0lz))  = E log  {p(y(t.)lY  .,0)}  + constant  (A. 7) 

i=l 

To  find  the  probability  distribution  of  y(t)  given  Y^_j^  and  0,  the  mean  value 
and  covariance  are  determined  first. 


E(y(t.)|Y._^,0)  A y(i/i-l) 


(A. 8) 


The  expected  value  or  the  mean  is  the  best  possible  estimate  of  measurements  at 
a point  given  the  measurements  until  the  previous  point. 


cov(y(t.)|Y._^.0)=  E{(y(t.)-y(i/i-l))(y(t.)-y(i/i-l))'^} 


A E{v(i)  v(i)M 


A B(i) 


(A. 9) 


v(i)  are  the  innovations  at  point  i and  B(i)  is  the  innovations  covariance. 

Since 

y(tp  - E(y(t.)|Y._j,0)  = v(i)  (A.  10) 

it  follows  that  y(h)  given  and  0 have  the  same  distribution  as  v(i). 

Kailath  [44  ] has  shown  that  as  the  sampling  rate  is  increased,  the  innovations 
v(i)  tend  towards  having  a Gaussian  density.  Assuming  a sufficiently  high 
sampling  rate,  the  distribution  of  v(i)  and,  therefore,  y(b)  given  and  0 

is  Gaussian,  i.e.  , 

expl-  y v(i)^  B ^(i)  v(i)i 

p(y(t  )|Y  ,0)  = L_2  ^ (A. 11) 

' ' ^ (27t)"’/2|B(i)|^/2 

log  {p(y(t )|Y^_j^,0)}  = - — v(i)  B ^(i)v(i)-i  log|B(i)|  + constant 


The  log-likelihood  function  of  Equation  (A. 7)  can  be  written  as 


log  (l^’Celz))  = - 7 £ {v'^(i)  B'^(i)  v(i)  + log|B(i)|) 

i=l 


An  estimate  of  the  unknown  parameters  is  obtained  by  maximizing  the  like- 
lihood function  or  the  log-likelihood  function  from  the  feasible  set  of  para- 


meter values. 
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(A.  14) 


9 = max  log(i2?(9|z.)) 
0e© 


= max  ■ T 2 {v^(i)B  ^(i)v(i)  + loglB(i)|} 
0e©  L i=l 


(A.  15) 


The  log -likelihood  function  depends  on  the  innovations  and  their  covar- 
iance. To  optimize  the  likelihood  function,  a way  must  be  found  for  determin- 
ing these  quantities.  Both  innovations  and  their  covariance  are  outputs  of  an 
extended  Kalman  filter.  This  required  Kalman  filter  is  developed  now. 


The  extended  Kalman  filter  is  conventionally  divided  into  two  parts . In 
the  first  part,  called  the  prediction  equations,  the  state  equations  and  state 
estimate  covariance  equations  are  propagated  in  time  from  one  measurement 
point  to  the  next.  In  the  second  part,  called  the  measurement  update  equations, 
the  measurements  and  associated  measurement  noise  covariances  are  used  to 
improve  state  estimates  . The  covariance  matrix  is  also  updated  at  this  point 
to  reflect  the  additional  information  obtained  from  the  measurements. 

Prediction  Equations 

The  state  prediction  is  done  using  the  equations  of  motion  (A. 2). 

Starting  at  time  with  current  estimate  x(i-l|i-l)  of  the  state  x(t^  and  the 
covariance  P(i-l|i-l),  the  following,  equations  are  used  to  find  the  predicted 
state  x(iji-l)  and  the  associated  covariance  P(iji-l);  see  Bryson  and  Ho  [ 45  ] . 

= f(x(tlt._^),uCt),0,t)  (A.  16) 


P(t|t._p  = F(t)  P(t|t._j)  + P(t|t._p  p'^ct)  + rqr'’^ 


t.  . < t < t. 
1-1  - - 1 


(A.  17) 
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The  nxn  matrix  F is  obtained  by  linearizing  f about  the  best  current 
estimate 


F(t)  = 


6f(x(t|t._p  ,u(t)  ,0,t) 

^x(t|t._^) 


(A.  18) 


Using  (A.  16)  to  (A.  18),  we  can  obtain 


x(t^lt._j^)  A x(ili-l) 
and  PCt.lt._p  A P(i]i-1) 


CA.19) 


Thereafter,  the  measurement  update  equations  are  used. 

Measurement  Update  Equations 

The  covariance  and  state  estimate  are  updated  using  the  measurements . 
The  necessary  relations  are  derived  by  Bryson  and  Ho  [45  ] and  are  presented 
here  without  proof.  The  innovation  and  its  covariance  are 


v(i)  = y(i)  - h(x(i|i-l),  u(b),  9,  t^ 
B(i)  = H(i)  P(i|i-1)  h’^Q)  + R 


(A. 20) 
(A. 21) 


where  H is  obtained  by  linearizing  h. 


H(i)  = 


dh(x  (i|  i-1)  ,u(t.  ,0,U) 
9x(ili-l) 


(A. 22) 


The  Kalman  gain  and  the  state  update  equations  are 


C 

* 

w ^ > 


K(i)  = P(iji-l)  B"^(i) 


x(ili)  = x(ili-l)  + K(i)  v(i) 


(A. 23) 
(A. 24) 
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w 


Finally,  P(i|i),  the  covariance  of  error  in  updated  state,  is  obtained  by 

P(i|i)  = (I-K(i)  H(i))  P(i|i-1)  (A. 25) 

One  can  now  return  to  the  time  update  or  prediction  equations . 

I 

j A.  3 OPTIMIZATION  PROCEDURE 

i Many  possible  numerical  procedures  can  be  used  for  this  optimization  prob 

; lem.  Modified  Newton-Raphson  [20]  or  Quasilinearization  [21]  have  been  found 

by  experience  to  give  quicker  convergence  than  most  procedures  like  the  con- 
jugate gradient  or  the  Davison  method.  The  modified  Newton-Raphson  is  a 
second  order  gradient  procedure  requiring  computation  of  first  and  second 

f 

order  partials  of  the  log -likelihood  function. 


Also 


The  last  three  terms  in  the  equation  for  second  partial  of  the  log-likelihood 
function  involve  second  partials  of  innovation  and  its  covariance.  Those  terms 
are  usually  dropped.  So  the  second  partial  is  approximated  by 


3^  log  (^(0|z)) 

3e.  36, 

J k 


^ |3v'^(i)  -1.  . 3v(i)  3v'^(i)  -L  3B(i)  -1.  . . . 

- s “^e; 36~® 

1=1  j k k j 


B-1(0  ^ B-l(i)  v(i) 
j k 


+ v^(i)B”^(i)  B'^(i)  ?|lilB“^(i)v(i) 

i k 


1 ^ L-l,.^  3B(i)  3B(i)l  I 

- ^Trj^B  (i)  -^B  (1)  -g^j  j 


(A. 28) 


The  gradients  of  innovation  and  its  covariance  for  parameter  0^  are  (Equations 
(A. 20)  and  (A. 21)) 


3v(i)  _ 3h| 


’ ^^lx=x(i|i-l) 


3x(i  i-1)  3h 
30  ~ W 

j = 1,2, ... ,p  : i = 1,2, 


(A. 29) 


Paii-l)H'^(i)  + H(i)  H^(i) 


3H"(i)  . 3R 

H(l)  P(l|l-1)  -go  + gg 

3 3 


(A. 30) 


j = 1,2,  ...  ,p 


i = 1.2 N 


Recursive  equations  can  be  obtained  for  gradients  of  the  predicted  state 
and  its  covariance.  This  is  done  in  stages  by  using  the  prediction  and 
measurement  update  equations  of  Kalman  filter.  Differentiating  (A.  16)  to  (A.  18) 
with  respect  to  0. 
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),u(t),t,0)  df(x(Wt.  ,),u(t),t,0) 


dx(t|t^_P 


9x  (0) 

ox(0|0)  _ o 


, ap(t|t.  ,)  .p...  ap(t|t,_ ,)  ap(t|t.  .)  „ 

ar  Vi>  * f 

J J J J 


prtit  '>  or^  + r E2.  + ro 


dP(0|0)  _ 


t.  1 < t < t 
1-1  — — 


ap(t)  ® f(x(t|t_p ,u(t) , 0,t) 
90j  3x 


j = 1,2, ...  ,p 


The  sensitivity  functions  are  updated  at  measurement  points  by  differentiating 
(A. 22)  to  (A. 25)  with  respect  to  0.. 


..  a'^h(x(ili-l)  ,u(tp,0,tp 


j 90j  ax(i|i-l) 


9KCO  _ H^(i)  B ^(i)  + P(i|i-1) 

90.  90. 

3 3 


-P(i|i-1)  h'^O)  B'^(i)  B“^(i 


(A. 36) 


0x(i|i) 

00. 

J 


0x(i|i-l) 

00. 

J 


0K(i) 
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J 


v(i)  + KCi) 


0v(i) 

00. 

J 


0P(i|i) 

00. 

J 


(I-K(i)  H(i)) 


0P(i  i-1) 

00. 

J 


0K(i) 

"5^ 

J 


H(i)  P(ili-l) 


- K(i)  P(i|i-1) 

j 


(A.  37) 


j = 1,2,3,  ...  ,p 


The  negative  of  the  matrix  of  second  partials  of  the  log-likelihood  func- 
tion is  called  the  information  matrix  M.  The  step  size  A0  for  parameter  esti- 
mates is  given  by 

ae  = ^ 

ou 

The  information  matrix  provides  a lower  bound  on  parameter  estimate  covar- 
iances , i . e . , 

E(0-0)  (O-O)"^  > m“^  (A. 39) 

V 

This  is  called  the  Cramer-Rao  lower  bound.  [46,47]  The  maximum  likelihood 
estimates  approach  this  bound  asymptotically. 

A.  4 UNEAR  SYSTEMS 

In  a linear  system,  the  functions  f and  h are  defined  as 


f(x,u,0,t)  AF(0,t)x  + G(0,t)u 
and  h(x,u,0,t)  A H(0,t)x  + D(6,t)u 


(A. 40) 
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In  this  case,  it  has  been  shown  that  if  the  model  is  correct,  the  innovations 
are  white  and  have  Gaussian  density  at  the  true  values  of  the  parameters . 
Therefore,  the  assumption  of  fast  sampling  rate  is  not  necessary. 

The  basic  algorithm  is  the  same.  However,  some  of  the  equations  can 
now  be  simplified.  The  equivalence  between  F(t)  and  H(i)  of  (A.  18)  and 
(A. 22)  and  F(0,t)  and  H(9,tp  of  (A. 40)  is  obvious.  Equations  (A.  16)  and 
(A. 20)  now  become 

^x(tjt._p  = F(0,t)  x(t|t._j^)  + G(0,t)  u(t)  (A. 41) 

and 

v(i)  = y(i)  - H(0,h)  x(i|i-l)  - D(0,t)  u(t.)  (A. 42) 

Equations  (A. 29)  and  (A. 31)  can  be  written  as 


= - H(0,t.)  gQ 
J 1 


0H(0,t.) 
dx(.ll-l) 1_ 


30. 

J 


3D(0,t.) 

~30; 

J 


and 


dT  — 30^ — 30^  ^ 30. 


dt  30. 
J 


j 


(^.43) 


(A.  44) 


All  other  equations  remain  the  same. 

There  is  considerable  reduction  in  computation  requirement  for  time- 
invariant  linear  system.  In  this  case,  matrices  F,  G,  H,  D,  P,  Q and  R and 


their  derivatives  with  respect  to  parameters  are  constant.  Gupta  [ 48  ] has  shown 
that  the  computation  of  state  sensitivities  can  be  reduced  to  many  fewer  equa- 
tions . Similar  reductions  are  possible  in  the  computation  of  covariance 
sensitivities . 

A. 5 TIME  INVARIANT  LINEAR  SYSTEMS  IN  STATISTICAL  STEADY  STATE 

In  many  aircraft  applications , the  Kalman  filter  is  in  steady  state  for  the 
duration  of  the  experiment.  This  occurs  when  the  Kalman  filter  is  in  opera- 
tion for  a sufficiently  long  time  and  the  process  and  measurement  noise  covari- 
ances do  not  change.  The  Kalman  gain  and  the  innovations  and  the  state 
covariances  approach  constant  values.  The  time  update  and  measurement  up 


date  equations  for  the  covariances  are 

|^P(t|t._p  = FP(tjt._^)'+  P(t|t._^)F'^  + (A. 45) 

K = Pali-DH"^  b"^  (A. 46) 

B = HP(i|i-l)H^  + R (A.  47) 

P(i|i)  = (I-KH)  P(i|i-1)  (A. 48) 

By  definition  of  the  steady  state 
P(i-l|i-l)  = P(i|i) 


Therefore,  from  (A. 45) 


P(i|i-U  = Pa-i|i-U  */‘ 


‘i-1 


(A.  49) 
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A (p(At)  (I-KH)  P(ili-l)  (p(At)  + Q' 


= cp(At)  CP(i|i-l)  - KBK  ) (p(At)  + Q' 


(A.  49) 


Using  (A. 46),  (A. 47)  and  (A. 49),  we  can  solve  for  P(i|i-1)  and  then  find 
K and  B.  Also,  it  can  be  shown  that 


3P  3P  T T 

II:  = a[  ^ - cpPA3Pcp^ 

J J 


(A. 50) 


==  (I-KH)  H^CHPH^  + R)‘^  + A. 

j j 


3B  _ 3H  p„T  . „ 3P  , „p  3H^  . 3R 
J J J J J 


(A. 51) 


(A.  52) 


where 


«> 

t 

* 

V ^ 

i*  ^ ^ 


1.  H 


Aj  = cp(I-KH) 

«-KH)P<p'^  * <p(I-KH)P  - <PK  II:  Pv"^  * !§-' 


J 


J 


J 


*3 ' Ir  If:  Ir  * if" 

j J J ) 


P A P(ili-l) 


3P  3K  3B 

Thus,  it  is  possible  to  solve  for  using  (A. 50)  and  then  find  ■ and 

j 3P  j ^ 

from  (A. 51)  and  (A. 52).  (A. 50)  is  a linear  equation  in  and  the 

j 

coefficient  of  the  unknown  matrix  does  not  depend  on  the  parameter  0..  Thus, 
the  sensitivity  of  state  covariance  matrix  can  be  determined  very  quickly  for 
all  parameters . Once  the  sensitivity  of  P , K and  B for  unknown  parameters 
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f 


is  determined,  only  state  sensitivity  equations  need  to  be  updated.  The  compu- 
tation of  state  sensitivity  functions  can  be  reduced  to  many  fewer  equations  as 
pointed  out  in  Section  A. 4. 

An  approximation  suggested  by  Mehra  [34  ] simplifies  the  problem  further. 
The  unknown  parameters  are  defined  to  include  elements  in  K and  B matrices 
instead  of  Q and  R.  Optimizing  the  log-likelihood  function  for  parameters  in 
B gives 


* 1 T 

B = i Z v(i)  v^(i) 

^ i=l 


(A. 54) 


The  gradient  of  the  log-likelihood  function  with  respect  to  other  unknown 
parameters  is 


i=l 


(A. 55) 


The  sensitivity  of  innovations  to  parameters  is  determined  using  the  following 
recursive  equations 

^x(t|t._j)  = Fx(t|t._p  + Gu(t) 


d 9x(t|t._^)  gp  9x(t|t._P  g^ 

dt  sw. — = m:  ^ “3^: — ^ m, 


j = 1,2,  ...  ,p  t 


v(i)  = y(i)  - Hx(i|i-1)  - Du(t^ 


x(i[i)  = xCiji-l)  + Kv(i) 
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Note  that 


=0  if  9.  is  not  an  element  of  K matrix 
00.  1 
J 

= I..,  , if  e.  A K..,  . 
j k J - J k 


(A. 59) 


where  is  a matrix  of  all  zeroes  except  a 1 at  the  j*,k*  position. 

This  approximation  simplifies  the  optimization  considerably.  However, 
this  usually  leads  to  an  overparameterized  model.  In  other  words,  once  K cind 
B are  determined , it  is  not  possible  to  find  any  corresponding  F , Q and  R 
which  have  the  desired  (known  a priori)  structure.  A good  estimate  of  ele- 
ments in  r,  Q and  R matrices  can  be  obtained  by  a least-squares  type  approach, 
The  fit  to  the  observed  data  i5  better  than  with  true  values  of  parameters 
but  the  parameter  estimatei'  do  not  have  minimum  variance.  This  approximation 
is  not  good  in  aircraft  application  where  the  structure  of  Q and  R is  known 
fairly  well,  but  is  excellent  where  there  are  many  process  noise  sources  and 
the  characteristics  of  both  Q and  R are  relatively  unknown  (e.g.,  economic 
systems)  . 

A. 6 MAXIMUM  LIKELIHOOD  WITH  NO  PROCESS  OR  MEASUREMENT  NOISE 


The  maximum  likelihood  method  can  be  simplified  when  either  process 
noise  or  measurement  noise  are  absent. 


9 
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No  Process  Noise 


»> 

t 


n' 


If  the  process  noise  is  zero  and  initial  states  are  known  perfectly,  i.e., 
w(t)  and  P(O)  are  zero,  the  covariance  of  the  error  in  the  predicted  state 
is  also  zero.  It  is  clear  from  (A. 23)  that  Kalman  gains  are  zero.  The  inno- 
vations is  the  output  error,  i.e.. 


v(i)  = y(i)  - h(x(tp  ,u(tp  ,0,t) 


(A. 60) 


and  the  innovation  covariance  is  (A. 21) 
B(i)  = R 

the  log-likelihood  function  is , 


log  0/f(6|z))  = - T 2 v’^(i)  R‘^v(i)  + log  |Rl 

i=l 

which  on  optimizing  for  unknown  parameters  in  R gives 


* 1 T 

R = i S vCi)  v^i) 

‘ i=l 


(A. 61) 


(A.  62) 


(A. 63) 


The  equality  in  (A. 63)  holds  only  for  those  elements  of  R which  are  not 
known  a priori.  For  instance,  even  if  R is  known  to  be  diagonal,  the  right 
hand  side  matrix  will  not  be  diagonal  in  general,  but  the  off-diagonal  terms 
should  be  ignored  before  they  are  equated  to  R.  Using  (A. 63)  in  (A. 62) 


logc^cejz)) 


1 N 


v^(i)  R"^v(i) 


+ constant 


(A. 64) 


The  optimizing  function  is  the  same  as  that  for  the  output  error  method  except 
that  the  measurement  noise  covariance  matrix  is  determined  using  (A. 63) 
and  is  used  as  the  weighting  matrix  in  the  criterion  function.  In  the  output 
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error  method,  the  measurement  noise  is  assumed  known  and  the  weighting 
function  is  arbitrary . 


(• 

¥ 

A- 

t 

i 

% 
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The  first  and  second  derivatives  of  the  log-likelihood  function  with 
respect  to  unknown  parameters  are 


^log(^(0|z))  = - Z v'^(i) 


i=l 


99. 

J 


(A. 65) 


log(  5^(01  z))  _ 

"90  96/ 

3 k 


J9v^(i)  A-1  3v(i) 

.r,)“55r“  ^ “§¥“ 

i=H  k ] 


i-1  3^v<l)  I 
) 


+ v"  (i)  R 


90.  96, 


(A. 66) 


The  terms  in  the  second  derivative  are  approximated  as 


9^  log(32’(0|z)) 

90.  99, 

3 k 


^ j 9v'^(i)  9v(i) 

1=1  1 k j 


(A.  67) 


No  Measurement  Noise 

If  all  states  are  measured  with  no  noise,  the  covariance  of  the  error  in 
state  estimates  is  zero  at  the  beginning  of  any  time  update, 


P(i-l|i-l)  = 0 


and  x(i-l|i-l)  = x(i-l) 


(A. 68) 


It  is  easy  to  show  in  this  case  that  for  fast  sampling  the  log-likelihood  function 
is  quadratic  in  the  difference  between  measured  values  of  x and  f(x,u,0,t). 

The  method  reduces  to  the  equation  error  method,  the  weight  W being  chosen 
as 
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Thus,  the  maximum  likelihood  method  and  equation-error  methods  are  equiva- 
lent except  for  the  technique  for  choosing  the  weighting  matrix. 


APPENDIX  B 

EQUATIONS  OF  MOTION  AND  LIST  OF  PARAMETERS 


Longitudinal  Case 


a 


state 

vector 


u 

q 

9 


L\J 


X * Fx  + Gu  -i-  Fw 
(angle-of-attack  perturbation) 
(longitudinal  velocity  perturbation) 
(pitch  rate  perturbation) 

(pitch  attitude  perturbation) 
(angle-of-attack  due  to  gusts) 


control 

vector 


(horizontal  elevator  deflection  from  trim) 


z 

u 

Z 

q 

-g/v  sin 

®o 

X 

u 

X 

q 

-g  cos 

0 

o 

F « 

\ 

M 

u 

M 

q 

0 

0 

0 

1 

0 

_ 0 

0 

0 

0 

’^5  ' 
e 

"0  “ 

e 

0 

G - 

e 

Fw  » 

0 

0 

0 

_ 0 

-"6- 
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whe  re 

g is  gravitational  acceleration, 

V is  steady-state,  total  velocity, 
9^  is  steady-state  pitch  angle 


measurements  y ■ Hx  + Du  + v 


Hx+Du 


q 

0 


-vfz  (a+a  )+Zr'5 

q <5^  e 


L"5. 


where 

i 

K = anele-of-attack  vane  scale  factor 
a ® 

I = longitudinal  distance  of  angle-of-attack  vane  forward 
of  c.g. 


Lateral  Case 


state 

vector 


(roll  rate) 

(yaw  rate) 

(sideslip  angle) 

(roll  angle) 

(sideslip  angle  due  to  gusts) 
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r 


control 

vector 


u = 


(aileron  deflection) 


(rudder  deflection) 
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* Not  identified 

[Ml  in  units  of  ft.,  sec.,  rad.) 

[Data  Supplied  by  NATC  to  SCI 

* Not  identified 

[All  in  unit-s  of-ft^  , rad.  ,-sec;) 

[Data  Supplied  by  NATC  TO  SCI  - July,  1974] 


* Not  identified 

[All  in  units  of  ft.,  rad.,  sec.) 

[Data  Supplied  by  NATC  TO  SCI  - April,  1974] 
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Table  3.3 
F-14  RUN2 

Altitude  30,000  ft. 
Speed  362  knots 
Sweepback  20° 


DERIVATIVE 

ST.ARTING  VALUE 
(LEAST  SQUARES 
'.AND  WIND  TUN^D 

1 

ESTIMATED 

RAW  DATA 

VALUE  imTH  la  BOUND) 
OUTLIERS  REMOVED 

z 

ct 

-0.38 

-0.347  (.033) 

-0.339  (.002"1) 

z 

-0.00066 

ft 

ft 

u 

z 

1.04 

ft 

A 

q 

0.0 

ft 

* 

X 

285.0 

* 

* 

a 

X 

-3.308 

ft 

ft 

u 

"q 

42.34 

* 

* 

M 

Ct 

-0.755 

-1.276  (.0789) 

-1.285  (.00649) 

M 

0.00416 

A 

it 

u 

-0.600 

-0.736  (.139) 

-0.7068  (.011) 

(jk) 

0.50 

* 

* 

c 

^<5e 

-0.20 

0.0378  (.218) 

0.119  (.020) 

X6e 

0.0 

A 

A 

^e 

-5.46 

-5.357  (.452) 

-5.338  (.033) 

Q 

0.0 

A 

* 

k 

1.0 

A 

* 

a 

kaila 

V 

0.0 

A 

-0.0754  (.0062) 

b 

0.0 

A 

0.000124 

a 

b 

0.0 

A 

-0.0934 

u 

b 

0.0 

A 

0.00056 

q 

^9 

0.0 

A 

0.000213 

0.0 

A 

-0.05 

a 

0.0075 

0.0428 

0.00159 

a 

0 

100.0 

52.0 

4.72 

u 

a 

0.007 

0.0424 

0.0038 

q 

0.0021 

0.015 

0.0877 

0.9 

8.1  . 

0.77 

az 

Table  3.4 


RUNl  (DOUBLET  INPUT) 

Altitude  30,000  ft. 
Speed  360  knots 
Sweepback  20® 


P.4RAVETER 

STARTING  VALUE 
(FROM  PREVIOUS  RUN) 

q FRCM  DATA 
(INCORREa) 

q NOT  USED 

q CORRECTED 
AND  USED 

-0.3394 

0.545 

-0.4268  (.0097) 

-0.4089  (.0059) 

L 

-0.00066 

* 

A 

A 

u 

Z 

1.04 

* 

A 

* 

q 

^0 

0.0 

* 

A 

ft 

X 

285.0 

* 

A 

A 

a 

X 

-3.3 

A 

A 

ft 

u 

X 

42.3 

A 

A 

it 

q 

M 

a 

M 

-1.28 

-4.59 

-2.002  (.018) 

-1.88  (.0078) 

0.0042 

* 

* 

A 

u 

-0.74 

-2.73 

-0.4907  (.024) 

-0.6726  (.0089) 

0.5 

A 

* 

* 

c 

^^e 

0.119 

0.639 

0.0856  (.014) 

-0.0147  (.0136) 

X6e 

0.0 

A 

M6e 

5.34 

0.416 

4.154  (.093) 

5.30  (.033) 

Q 

0.0 

A 

* 

A 

k 

1.0 

A 

A 

0.881  (.0089) 

a 

kaila 

V 

0.0 

-0.08 

-0.0594  (.0066) 

A 

b 

0.0 

-0.0328 

0.00609 

-0.00071 

a 

b 

0.0 

-2.06 

0.344 

-0.147 

u 

b 

0.0 

-0.0104 

0.00164 

-0.000877 

q 

-0.00505 

baj 

0.0 

-0.0378 

0.0112 

0.0 

-0.81 

0.922 

-0.528 

0 

0.0075 

0.034 

0.0072 

0.00262 

a 

a 

100.0 

2.7 

4.05 

1.35 

u 

a 

0.007 

0.023 

0.054 

0.0034 

q 

0.0100 

0.015 

0.063 

0.0131 

0.9 

10.0 

2.6 

2.1 

Table  3.5 


LATERAL  MOTIONS  (RUNS)* 

Altitude  30,000  ft. 
Speed  360  knots 
Sweepback  20° 


1 

STARTING  value' 
(LEAST  SQUARES) 

MAXIMUM  LIKELIHOOD  ESTIMATE  (WITH  la  BOUND) 

PARAMETER 

NO  PROCESS  NOISE 

PROCESS  ?«3ISE  

i 

Lp 

-1.975 

-2.004  (.0235) 

-1.948  (.0205) 

Lr 

4.877 

2.686  (.0645) 

2.369  (.0541) 

‘•3 

-8.29 

-8.500  (.0982) 

-6.952  (.0846) 

N 

-0.004 

* 

A 

P 

-0.62 

* 

A 

1.96 

* 

A 

-0.05 

-0.0649  (.00190) 

-0.0819  (.00121) 

a 

0.0 

it 

A 

0 

9 

0.0 

* 

A 

0 

L6a 

10.29 

* 

A 

^5a 

0.24 

A 

it 

^^a 

0.0 

* 

A 

-2.828 

-3.337  (.0865) 

-1.446  (.0661) 

1.55 

A 

A 

0.42 

-0.015  (.0018) 

-0.0150  (.00121) 

1.0 

» 

A 

Vs^ 

0.03 

A 

A 

0) 

0.5 

A 

* 

c 

Q 

0.0 

A 

0.000547  (.000026) 

b 

-0.000209 

-0.00264 

P 

- 

-0.000126 

0.000318 

»8 

tay 

- 

-0.ofif.a590 

0.000636 

- 

-0.00230 

-0.00516 

- 

0.00568 

-0.0629 

a 

0.005 

0.0405 

0.0246 

p 

a 

0.005 

0.0196 

0.0069 

^av 

/ 

0.005 

0.0190 

0.00487 

0.005 

0.0192 

0.0245 

0.5 

1.57 

0.85 

